## onefile=FALSE
## pdf(file = ifelse(onefile, "Rplots.pdf", "Rplot%03d.pdf"),
##         encoding="ISOLatin1",
##         width = 10, height = 8, onefile = TRUE, family = "Helvetica",
##         title = "R Graphics Output", fonts = NULL, version = "1.1")

##
## setting up plotting space
##
par(mfrow = c(1, 2))

##
## reading data and setting parameters for the cascade (left) panel
##
data_cascade <- read.csv(file="../figure2_cascade.csv",header=FALSE)
title_cascade <- "Cascade routine is ON"

##
## reading data and setting parameters for the NOcascade (right) panel
##

data_nocascade <- read.csv(file="../figure2_nocascade.csv",header=FALSE)
title_nocascade <- "Cascade routine is OFF"


##
## left panel plotting section
##
#
# margins
par(mar = c(5, 5, 3, 0))
#
# data dimension
rows_num <- dim(data_cascade)[1]
#
# min and max axis values for X and Y
minX.pow <- 3
maxX.pow <- 5
minY.pow <- -3
maxY.pow <- 0
# draw an empty plot frame
plot(NULL, xlim=c(10^minX.pow, 10^maxX.pow), ylim=c(10^minY.pow, 10^maxY.pow),
     log = "xy", xaxt = "n",  yaxt = "n", type = "l", xlab=expression(paste(n[0](cm^{-3}))), 
     ylab=list("Fractional abundance", cex=1.2), col=1)
# add the title
title(main=list(title_cascade, cex=0.9))
# draw all the abundances lines
for(i in 0:rows_num) {
 dataToPlot <- t( rbind(c(1E3,1E4,1E5,1E6), data_cascade[i,2:5]) )
 lines(dataToPlot, col=i, lty=i, lwd=1.5)
 points(dataToPlot, pch=i, cex=1.5, col=i, lty=i)
}
# set-up labels
labels.text=t(data_cascade)[1,]
labels.col=seq(1,rows_num,1)
# set the axis major tick marks
axisX.at <- 10 ^ c(minX.pow:maxX.pow)
axisY.at <- 10 ^ c(minY.pow:maxY.pow)
# draw the major tick marks and label them using plotmath 
axis(1, at = axisX.at, tcl = -1, labels = parse(text = paste("10^", minX.pow:maxX.pow, sep = "")))
axis(2, at = axisY.at, tcl = -1, labels = parse(text = paste("10^", minY.pow:maxY.pow, sep = "")))
# now do the minor ticks, at 1/10 of each power of 10 interval 
axis(1, at = 1:10 * rep(axisX.at[-1] / 10, each = 10), tcl = -0.5, labels = FALSE)
axis(2, at = 1:10 * rep(axisY.at[-1] / 10, each = 10), tcl = -0.5, labels = FALSE)
#add the legend
legend(1.1E3,0.22,labels.text, col=labels.col, 
		text.col = "green4", lty = labels.col, 
		pch = labels.col, ncol=1, merge = TRUE, cex=0.8)

##
## right panel plotting section
##
#
# margins
par(mar = c(5, 0, 3, 3))
#
# data dimension
rows_num <- dim(data_nocascade)[1]
#
# min and max axis values for X and Y
minX.pow <- 3
maxX.pow <- 5
minY.pow <- -3
maxY.pow <- 0
# draw an empty plot frame
plot(NULL, xlim=c(10^minX.pow, 10^maxX.pow), ylim=c(10^minY.pow, 10^maxY.pow),
     log = "xy", xaxt = "n",  yaxt = "n", type = "l", xlab=expression(paste(n[0](cm^{-3}))), col=1)
# add the title
title(main=list(title_nocascade, cex=0.9))
# draw all the abundances lines
for(i in 0:rows_num) {
 dataToPlot <- t( rbind(c(1E3,1E4,1E5,1E6), data_nocascade[i,2:5]) )
 lines(dataToPlot, col=i, lty=i, lwd=1.5)
 points(dataToPlot, pch=i, cex=1.5, col=i, lty=i)
}
# set-up labels
labels.text=t(data_nocascade)[1,]
labels.col=seq(1,rows_num,1)
# set the axis major tick marks
axisX.at <- 10 ^ c(minX.pow:maxX.pow)
axisY.at <- 10 ^ c(minY.pow:maxY.pow)
# draw the major tick marks and label them using plotmath 
axis(1, at = axisX.at, tcl = -1, labels = parse(text = paste("10^", minX.pow:maxX.pow, sep = "")))
axis(4, at = axisY.at, tcl = -1, labels = parse(text = paste("10^", minY.pow:maxY.pow, sep = "")))
axis(2, at = axisY.at, tcl = -0.5, labels=FALSE)
axis(2, at = axisY.at, tcl = 0.5, labels=FALSE)
# now do the minor ticks, at 1/10 of each power of 10 interval 
axis(1, at = 1:10 * rep(axisX.at[-1] / 10, each = 10), tcl = -0.5, labels = FALSE)
axis(4, at = 1:10 * rep(axisY.at[-1] / 10, each = 10), tcl = -0.5, labels = FALSE)
axis(2, at = 1:10 * rep(axisY.at[-1] / 10, each = 10), tcl = -0.2, labels=FALSE)
axis(2, at = 1:10 * rep(axisY.at[-1] / 10, each = 10), tcl = 0.2, labels=FALSE)
#add the legend
legend(1.1E3,0.22,labels.text, col=labels.col, 
		text.col = "green4", lty = labels.col, 
		pch = labels.col, ncol=1, merge = TRUE, cex=0.8)

##
## tear down plotting
## dev.off()